rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)

#This is to obtain:
#-figure_a11a
#-figure_a11b
#-figure_a11c
#-figure_a11d


library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
library("corrplot")
library("Hmisc")
library("glue")
library("ggplot2")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library("plm")
library("multiwayvcov")
library("lmtest")
library("psych")
library("ggrepel")
library("lemon")
library("sandwich")
library("spdep")
library("xtable")
library('gsynth')
library("pracma")
library("conleyreg")
library("vars")
library("rlist")
library("pgirmess")
library("geojsonsf")

#Step1. Turning some variables to numeric
data_cov_mob = geojson_sf("./data/data_cov_mob.geojson")
data_cov_mob$pct_manufacturing<-as.numeric(data_cov_mob$pct_manufacturing)
data_nocovid_nomob<-subset(data_cov_mob, week=="00")

data_cov_mob$kul<-as.numeric(data_cov_mob$kul)
data_cov_mob$nat<-as.numeric(data_cov_mob$nat)
data_cov_mob$frei<-as.numeric(data_cov_mob$frei)
data_cov_mob$soz<-as.numeric(data_cov_mob$soz)
data_cov_mob$pol<-as.numeric(data_cov_mob$pol)
data_cov_mob$inter<-as.numeric(data_cov_mob$inter)

data_nocovid_nomob$kul<-as.numeric(data_nocovid_nomob$kul)
data_nocovid_nomob$nat<-as.numeric(data_nocovid_nomob$nat)
data_nocovid_nomob$frei<-as.numeric(data_nocovid_nomob$frei)
data_nocovid_nomob$soz<-as.numeric(data_nocovid_nomob$soz)
data_nocovid_nomob$pol<-as.numeric(data_nocovid_nomob$pol)
data_nocovid_nomob$inter<-as.numeric(data_nocovid_nomob$inter)

data_cov_mob$gspo<-as.numeric(data_cov_mob$gspo)
data_cov_mob$gkul<-as.numeric(data_cov_mob$gkul)
data_cov_mob$gnat<-as.numeric(data_cov_mob$gnat)

data_cov_mob$gsoz<-as.numeric(data_cov_mob$gsoz)
data_cov_mob$gpol<-as.numeric(data_cov_mob$gpol)
data_cov_mob$ginter<-as.numeric(data_cov_mob$ginter)

#Bridging Networks
data_cov_mob$bridging_tot_pop<-data_cov_mob$bridging_total/data_cov_mob$pop_total*1000
mean_brid<-summary(data_cov_mob$bridging_tot_pop)[4]
median_brid<-summary(data_cov_mob$bridging_tot_pop)[3]
upper_quart_brid<-summary(data_cov_mob$bridging_tot_pop)[5]
data_cov_mob$bridging_tot_pop_dummy_mean<-ifelse(data_cov_mob$bridging_tot_pop>mean_brid, 1, 0)
data_cov_mob$bridging_tot_pop_dummy_med<-ifelse(data_cov_mob$bridging_tot_pop>median_brid, 1, 0)
data_cov_mob$bridging_tot_pop_dummy_quar<-ifelse(data_cov_mob$bridging_tot_pop>upper_quart_brid, 1, 0)

data<-data_cov_mob
data$week_new<-as.numeric(data$week)+1
data$week_new<-paste0("0", data$week_new)  
data$week_new2<-ifelse(nchar(data$week_new)==3, str_remove(data$week_new, "^0+"), data$week_new)
data$week<-data$week_new2
data<-subset(data, select = -c(week_new, week_new2))
data$post_treat<-ifelse(as.numeric(data$week)>10, 1,0)
data$post_treat_bridging_tot_pop_dummy_mean<-data$post_treat*data$bridging_tot_pop_dummy_mean
data$post_treat_bridging_tot_pop_dummy_med<-data$post_treat*data$bridging_tot_pop_dummy_med
data$post_treat_bridging_tot_pop_dummy_quar<-data$post_treat*data$bridging_tot_pop_dummy_quar
data$post_treat_bridging_tot_pop_dummy_cont<-data$post_treat*data$bridging_tot_pop

###########
#All clubs#
###########

#Create dummies
library('fastDummies')
# Create dummy variables:
dataf <- dummy_cols(data, select_columns = 'AGS')
dataw <- dummy_cols(dataf, select_columns = 'week')
datas <- dummy_cols(dataw, select_columns = 'SN_L')

library(dplyr)
new<-dataw %>% dplyr::select(starts_with("AGS_"))
new2<-dataw %>% dplyr::select(starts_with("week_"))
new3<-datas %>% dplyr::select(starts_with("SN_L"))

list_dummies<-names(new)
#I am removing the counties which get dropped in Stata

#Re-enable this later if necessary
list_dummies<-list_dummies[list_dummies!= "AGS_0"]
list_dummies<-list_dummies[list_dummies!= "AGS_5316"]
list_dummies<-list_dummies[list_dummies!= "AGS_8121"]
list_dummies<-list_dummies[list_dummies!= "AGS_8221"]
list_dummies<-list_dummies[list_dummies!= "AGS_8231"]
list_dummies<-list_dummies[list_dummies!= "AGS_8311"]
list_dummies<-list_dummies[list_dummies!= "AGS_8421"]
list_dummies<-list_dummies[list_dummies!= "AGS_9461"]
list_dummies<-list_dummies[list_dummies!= "AGS_13004"]
list_dummies<-list_dummies[list_dummies!= "AGS_16052"]
list_dummies<-list_dummies[list_dummies!= "AGS_16053"]
list_dummies<-list_dummies[list_dummies!= "AGS_16054"]
list_dummies<-list_dummies[list_dummies!= "AGS_16055"]
list_dummies<-list_dummies[list_dummies!= "AGS_16061"]
list_dummies<-list_dummies[list_dummies!= "AGS_16077"]

list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]
#Re-enable this later if necessary
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_09"]
list_week_dummies2<-list_week_dummies2[list_week_dummies2!= "week_52"]


list_land_dummies<-names(new3)
#Re-enable this later if necessary
list_land_dummies_full<-list_land_dummies[list_land_dummies!= "SN_L"]
list_land_dummies<-list_land_dummies_full[list_land_dummies_full!= "SN_L_01"]

#Creating the interaction terms for bonding
for (i in list_land_dummies){
  for (j in list_week_dummies2)
    datas[[paste("inter_", i, "_", j, sep="")]]<-datas[[i]]*datas[[j]]
}

new<-datas %>% dplyr::select(starts_with("inter_"))
list_land_dummies_star<-names(new)

#Event study
#Bridging Interaction
new2<-dataw %>% dplyr::select(starts_with("week_"))
list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]

#Creating the interaction terms for bridging
for (i in list_week_dummies){
  datas[[paste("bri_mean_", i, sep="")]]<-datas[[i]]*datas$bridging_tot_pop_dummy_mean
  datas[[paste("bri_med_", i, sep="")]]<-datas[[i]]*datas$bridging_tot_pop_dummy_med
  datas[[paste("bri_quar_", i, sep="")]]<-datas[[i]]*datas$bridging_tot_pop_dummy_quar
  datas[[paste("bri_cont_", i, sep="")]]<-datas[[i]]*datas$bridging_tot_pop
  }

#Bri
inter_terms_bri<-datas %>% dplyr::select(starts_with("bri_mean_week"))
inter_terms_bri<-names(inter_terms_bri)
inter_terms_bri<-inter_terms_bri[inter_terms_bri != "bri_mean_week_10"]

inter_terms_bri<-datas %>% dplyr::select(starts_with("bri_med_week"))
inter_terms_bri<-names(inter_terms_bri)
inter_terms_bri<-inter_terms_bri[inter_terms_bri != "bri_med_week_10"]

inter_terms_bri<-datas %>% dplyr::select(starts_with("bri_quar_week"))
inter_terms_bri<-names(inter_terms_bri)
inter_terms_bri<-inter_terms_bri[inter_terms_bri != "bri_quar_week_10"]

inter_terms_bri<-datas %>% dplyr::select(starts_with("bri_cont_week"))
inter_terms_bri<-names(inter_terms_bri)
inter_terms_bri<-inter_terms_bri[inter_terms_bri != "bri_cont_week_10"]

######
#MEAN#
######

inter_terms_bri<-datas %>% dplyr::select(starts_with("bri_mean_week"))
inter_terms_bri<-names(inter_terms_bri)
inter_terms_bri<-inter_terms_bri[inter_terms_bri != "bri_mean_week_10"]

#Deleting variables which are collinear
list_land_dummies_full2<-list_land_dummies_full[list_land_dummies_full!= "SN_L_02"]
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_53"]
for_event_bri2<-list.append(inter_terms_bri, list_land_dummies_full2, list_week_dummies2)

#Deleting missing observations
datas$time<-as.numeric(datas$week)
case_bri<-subset(datas, select = list.append(for_event_bri2, "bridging_total", "mean_mobil", "POINT_X", "POINT_Y",
                                         "time", "AGS", "east_ger"))
case_bri<-case_bri[complete.cases(case_bri), ]

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(case_bri$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)

#Merging with spatial dataframe
nuts<-subset(data_cov_mob, week=="00")
nuts2<-subset(nuts, select = c("AGS"))
nuts2$AGS<-as.numeric(as.character(nuts2$AGS))
case_bri_oneweek<-subset(case_bri, week_11==1)

#Bridging
#Merging with spatial dataframe
nuts$AGS<-as.numeric(nuts$AGS)
nuts_merged<-left_join(nuts2, case_bri_oneweek, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))

#Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

#Getting coordinates
coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)

#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)
dist_cut_bri<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

#Bri formula
time_trend_bri<-as.formula(paste("mean_mobil~",
                                 paste(for_event_bri2, collapse= "+")))

#BRIDGING: Creating distance matrix - accelerates the conely regression
dm_bri <- dist_mat(case_bri, 
                   lat = "POINT_Y", lon = "POINT_X", 
                   unit = "AGS", time = "time")

df_bri<-conleyreg(time_trend_bri, data=case_bri, 
                  dist_cut_bri, unit = "AGS", time = "time", 
                  dist_mat=dm_bri, lag_cutoff=temp_lag)
df_bri<-broom::tidy(df_bri)
df_bri$term[grepl( "bri_mean_week_" , df_bri$term)]
df_bri<-add_row(df_bri,term="bri_mean_week_10",estimate=0)
df_bri<-subset(df_bri, grepl( "bri_mean_week_" , df_bri$term))
df_bri<-df_bri[order(df_bri$term),]

x<-seq(-10, 42, by = 1)
df_bri$observation<-NA
df_bri$observation <- x

df_bri$observation<-as.character(df_bri$observation)
selection<-df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]
df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]<-paste("+", selection, sep = "")
df_bri$t<-paste("Network x \n", "(t", df_bri$observation, ")", sep = "")
df_bri$t[df_bri$t=="Network x \n(t+0)"]<-"Network x t"

df_bri$week<-NA
df_bri$week<-readr::parse_number(df_bri$term)
df_bri$model<-"Bridging"


###########################
#Zoom Event Study Bridging#
###########################

df_bri_small<-subset(df_bri, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_bri_small[as.numeric(df_bri_small$week) %in% irislabs1,]$t)

coef_zoom_bri<-ggplot(df_bri_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_bri, file = "./paper/graphs/figure_a11a.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)

#########
#MEDIAN#
########

inter_terms_bri<-datas %>% dplyr::select(starts_with("bri_med_week"))
inter_terms_bri<-names(inter_terms_bri)
inter_terms_bri<-inter_terms_bri[inter_terms_bri != "bri_med_week_10"]

#Deleting variables which are collinear
list_land_dummies_full2<-list_land_dummies_full[list_land_dummies_full!= "SN_L_02"]
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_53"]
for_event_bri2<-list.append(inter_terms_bri, list_land_dummies_full2, list_week_dummies2)

#Deleting missing observations
datas$time<-as.numeric(datas$week)
case_bri<-subset(datas, select = list.append(for_event_bri2, "bridging_total", "mean_mobil", "POINT_X", "POINT_Y",
                                             "time", "AGS", "east_ger"))
case_bri<-case_bri[complete.cases(case_bri), ]

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(case_bri$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)

#Merging with spatial dataframe
nuts<-subset(data_cov_mob, week=="00")
nuts2<-subset(nuts, select = c("AGS"))
nuts2$AGS<-as.numeric(as.character(nuts2$AGS))
case_bri_oneweek<-subset(case_bri, week_11==1)

#Merging with spatial dataframe
nuts$AGS<-as.numeric(nuts$AGS)
nuts_merged<-left_join(nuts2, case_bri_oneweek, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))

#Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)

#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)
dist_cut_bri<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

#Bri formula
time_trend_bri<-as.formula(paste("mean_mobil~",
                                 paste(for_event_bri2, collapse= "+")))

#BRIDGING: Creating distance matrix - accelerates the conely regression
dm_bri <- dist_mat(case_bri, 
                   lat = "POINT_Y", lon = "POINT_X", 
                   unit = "AGS", time = "time")

df_bri<-conleyreg(time_trend_bri, data=case_bri, 
                  dist_cut_bri, unit = "AGS", time = "time", 
                  dist_mat=dm_bri, lag_cutoff=temp_lag)
df_bri<-broom::tidy(df_bri)
df_bri$term[grepl( "bri_med_week_" , df_bri$term)]
df_bri<-add_row(df_bri,term="bri_med_week_10",estimate=0)
df_bri<-subset(df_bri, grepl( "bri_med_week_" , df_bri$term))
df_bri<-df_bri[order(df_bri$term),]

x<-seq(-10, 42, by = 1)
df_bri$observation<-NA
df_bri$observation <- x

df_bri$observation<-as.character(df_bri$observation)
selection<-df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]
df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]<-paste("+", selection, sep = "")
df_bri$t<-paste("Network x \n", "(t", df_bri$observation, ")", sep = "")
df_bri$t[df_bri$t=="Network x \n(t+0)"]<-"Network x t"

df_bri$week<-NA
df_bri$week<-readr::parse_number(df_bri$term)
df_bri$model<-"Bridging"

###########################
#Zoom Event Study Bridging#
###########################

df_bri_small<-subset(df_bri, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_bri_small[as.numeric(df_bri_small$week) %in% irislabs1,]$t)

coef_zoom_bri<-ggplot(df_bri_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_bri, file = "./paper/graphs/figure_a11b.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)

##########
#QUARTILE#
##########

inter_terms_bri<-datas %>% dplyr::select(starts_with("bri_quar_week"))
inter_terms_bri<-names(inter_terms_bri)
inter_terms_bri<-inter_terms_bri[inter_terms_bri != "bri_quar_week_10"]

#Deleting variables which are collinear
list_land_dummies_full2<-list_land_dummies_full[list_land_dummies_full!= "SN_L_02"]
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_53"]
for_event_bri2<-list.append(inter_terms_bri, list_land_dummies_full2, list_week_dummies2)

#Deleting missing observations
datas$time<-as.numeric(datas$week)
case_bri<-subset(datas, select = list.append(for_event_bri2, "bridging_total", "mean_mobil", "POINT_X", "POINT_Y",
                                             "time", "AGS", "east_ger"))
case_bri<-case_bri[complete.cases(case_bri), ]
summary(case_bri$bri_mean_week_30)

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(case_bri$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)

#Merging with spatial dataframe
nuts<-subset(data_cov_mob, week=="00")
#nuts<-subset(data, select = c("AGS"))
nuts2<-subset(nuts, select = c("AGS"))
nuts2$AGS<-as.numeric(as.character(nuts2$AGS))
case_bri_oneweek<-subset(case_bri, week_11==1)

#Merging with spatial dataframe
nuts$AGS<-as.numeric(nuts$AGS)
nuts_merged<-left_join(nuts2, case_bri_oneweek, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))


#Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)

#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)
dist_cut_bri<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

#Bri formula
time_trend_bri<-as.formula(paste("mean_mobil~",
                                 paste(for_event_bri2, collapse= "+")))

#BRIDGING: Creating distance matrix - accelerates the conely regression
dm_bri <- dist_mat(case_bri, 
                   lat = "POINT_Y", lon = "POINT_X", 
                   unit = "AGS", time = "time")

df_bri<-conleyreg(time_trend_bri, data=case_bri, 
                  dist_cut_bri, unit = "AGS", time = "time", 
                  dist_mat=dm_bri, lag_cutoff=temp_lag)
df_bri<-broom::tidy(df_bri)
df_bri$term[grepl( "bri_quar_week_" , df_bri$term)]
df_bri<-add_row(df_bri,term="bri_quar_week_10",estimate=0)
df_bri<-subset(df_bri, grepl( "bri_quar_week_" , df_bri$term))
df_bri<-df_bri[order(df_bri$term),]

x<-seq(-10, 42, by = 1)
df_bri$observation<-NA
df_bri$observation <- x

df_bri$observation<-as.character(df_bri$observation)
selection<-df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]
df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]<-paste("+", selection, sep = "")
df_bri$t<-paste("Network x \n", "(t", df_bri$observation, ")", sep = "")
df_bri$t[df_bri$t=="Network x \n(t+0)"]<-"Network x t"

df_bri$week<-NA
df_bri$week<-readr::parse_number(df_bri$term)
df_bri$model<-"Bridging"


###########################
#Zoom Event Study Bridging#
###########################

df_bri_small<-subset(df_bri, week>=(5) & week<=(16))

irislabs1<-seq(5,16,1)
irislabs2<-unique(df_bri_small[as.numeric(df_bri_small$week) %in% irislabs1,]$t)

coef_zoom_bri<-ggplot(df_bri_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_bri, file = "./paper/graphs/figure_a11c.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)

############
#Continuous#
############

inter_terms_bri<-datas %>% dplyr::select(starts_with("bri_cont_week"))
inter_terms_bri<-names(inter_terms_bri)
inter_terms_bri<-inter_terms_bri[inter_terms_bri != "bri_cont_week_10"]

#Deleting variables which are collinear
list_land_dummies_full2<-list_land_dummies_full[list_land_dummies_full!= "SN_L_02"]
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_53"]
for_event_bri2<-list.append(inter_terms_bri, list_land_dummies_full2, list_week_dummies2)

#Deleting missing observations
datas$time<-as.numeric(datas$week)
case_bri<-subset(datas, select = list.append(for_event_bri2, "bridging_total", "mean_mobil", "POINT_X", "POINT_Y",
                                             "time", "AGS", "east_ger"))
case_bri<-case_bri[complete.cases(case_bri), ]

#Deciding on the temporal lag - value with the lowest AIC
x<-VARselect(case_bri$mean_mobil)
min(x$criteria["AIC(n)",])
temp_lag<-names(x$criteria["AIC(n)",])[x$criteria["AIC(n)",]==min(x$criteria["AIC(n)",])]
temp_lag<-as.numeric(temp_lag)

#Merging with spatial dataframe
nuts<-subset(data_cov_mob, week=="00")
nuts2<-subset(nuts, select = c("AGS"))
nuts2$AGS<-as.numeric(as.character(nuts2$AGS))
case_bri_oneweek<-subset(case_bri, week_11==1)

#Bridging
nuts_merged<-left_join(nuts2, case_bri_oneweek, by = c("AGS" = "AGS"))
nuts_merged2<-subset(nuts_merged, !is.na(nuts_merged$bridging_total))

#Redefinining projection to calculate dist in meters
nuts_merged3 <- st_transform(nuts_merged2, crs = "+proj=utm +zone=1 +units=cm")
nuts_merged3$POINT_X_m<-st_coordinates(st_centroid(nuts_merged3))[,1]
nuts_merged3$POINT_Y_m<-st_coordinates(st_centroid(nuts_merged3))[,2]

coo<-subset(nuts_merged3, select = c("POINT_X_m", "POINT_Y_m"))
coo<-st_drop_geometry(coo)


#Calculating the correlogram
pgi_cor <- correlog(coords=coo, z=nuts_merged3$bridging_total, method="Moran", 
                    nbclass=10)
pgi_cor_df<-data.frame(pgi_cor)
dist_cut_bri<-round(pgi_cor_df$dist.class[abs(pgi_cor_df$coef)== min(abs(pgi_cor_df$coef))]/100000,0)

#Bri formula
time_trend_bri<-as.formula(paste("mean_mobil~",
                                 paste(for_event_bri2, collapse= "+")))

#BRIDGING: Creating distance matrix - accelerates the conely regression
dm_bri <- dist_mat(case_bri, 
                   lat = "POINT_Y", lon = "POINT_X", 
                   unit = "AGS", time = "time")

df_bri<-conleyreg(time_trend_bri, data=case_bri, 
                  dist_cut_bri, unit = "AGS", time = "time", 
                  dist_mat=dm_bri, lag_cutoff=temp_lag)
df_bri<-broom::tidy(df_bri)
df_bri$term[grepl( "bri_cont_week_" , df_bri$term)]
df_bri<-add_row(df_bri,term="bri_cont_week_10",estimate=0)
df_bri<-subset(df_bri, grepl( "bri_cont_week_" , df_bri$term))
df_bri<-df_bri[order(df_bri$term),]

x<-seq(-10, 42, by = 1)
df_bri$observation<-NA
df_bri$observation <- x

df_bri$observation<-as.character(df_bri$observation)
selection<-df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]
df_bri$observation[!stringr::str_detect(df_bri$observation, '-')]<-paste("+", selection, sep = "")
df_bri$t<-paste("Network x \n", "(t", df_bri$observation, ")", sep = "")
df_bri$t[df_bri$t=="Network x \n(t+0)"]<-"Network x t"

df_bri$week<-NA
df_bri$week<-readr::parse_number(df_bri$term)
df_bri$model<-"Bridging"


###########################
#Zoom Event Study Bridging#
###########################

df_bri_small<-subset(df_bri, week>=(5) & week<=(16))

irislabs1<-seq(1,53,1)
irislabs2<-unique(df_bri[as.numeric(df_bri$week) %in% irislabs1,]$t)
length(irislabs2)

coef_zoom_bri<-ggplot(df_bri_small, aes(x=week, y=estimate))+
  geom_point(size = 2,
             position=position_dodge(width=0.4))+
  geom_errorbar(aes(ymin = estimate-1.96*std.error, ymax = estimate+1.96*std.error), 
                size = 1, width = 0,
                position=position_dodge(width=0.4))+
  scale_x_continuous(breaks = irislabs1,
                     labels = irislabs1,
                     sec.axis = sec_axis(~.,
                                         breaks =irislabs1,
                                         labels = irislabs2))+
  geom_hline(yintercept = 0) +
  geom_vline(xintercept=10,color="red")+
  #scale_x_continuous(name = "Bandwidth (km)", breaks=seq(5,85,5)) +
  labs(x = "Week", y = "Point Estimate")+
  theme_bw() +
  theme(axis.text.x = element_text(size=12),
        axis.text.x.top  = element_text(size=10, angle = 45, margin=margin(5,5,10,5,"pt")),
        axis.text.y = element_text(size=12),
        axis.title=element_text(size=12),
        plot.title = element_text(hjust = 0.5),
        legend.position = c(1, 0),
        #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
        #and c(1,1) corresponds to the "top right" position.
        legend.box.background = element_rect(fill='white'),
        legend.background = element_blank(),
        legend.text=element_text(size=12))
#plot_ex_coef2<-reposition_legend(plot_ex_coef, 'bottom left')
ggsave(coef_zoom_bri, file = "./paper/graphs/figure_a11d.jpg", 
       height = 15, width = 20, 
       units = "cm", dpi = 200)


